arXiv:1506.02583vl [math.OC] 8 Jun 2015 


Preconditioned Continuation Model Predictive Control 


Andrew Knyazev* Yuta Fujii^ Alexander Malyshev^ 


Abstract 

Model predictive control (MFC) anticipates future events to 
take appropriate control actions. Nonlinear MFC (NMFC) 
describes systems with nonlinear models and/or constraints. 
A Continuation/GMRES Method for NMFC, suggested by 
T. Ohtsuka in 2004, uses the GMRES iterative algorithm to 
solve a forward difference approximation Ax — b of the Con¬ 
tinuation NMFC (CNMFC) equations on every time step. 
The coefficient matrix A of the linear system is often ill- 
conditioned, resulting in poor GMRES convergence, slowing 
down the on-line computation of the control by CNMFC, 
and reducing control quality. We adopt CNMFC for chal¬ 
lenging minimum-time problems, and improve performance 
by introducing efficient preconditioning, utilizing parallel 
computing, and substituting MINRES for GMRES. 

1 Introduction 

Model predictive control (MFC) is used in many appli¬ 
cations to control complex dynamical systems. Exam¬ 
ples of such systems include production lines, car en¬ 
gines, robots, other numerically controlled machining, 
and power generators. The MFC is based on optimiza¬ 
tion of the operation of the system over a future finite 
time-horizon, subject to constraints, and implementing 
the control only over the current time step. 

Model predictive controllers rely on dynamic mod¬ 
els of the process, most often linear empirical models, in 
which case the MFC is linear. Nonlinear MFC (NMFC), 
which describes systems with nonlinear models and con¬ 
straints, is often more realistic, compared to the linear 
MFC, but computationally more difficult. Similar to 
the linear MFC, the NMFC requires solving optimal 
control problems on a finite prediction horizon, gener¬ 
ally not convex, which poses computational challenges. 
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Numerical solution of the NMFC optimal control 
problems may be based on Newton-type optimization 
schemes. Exact Newton-type optimization schemes re¬ 
quire an analytic expression of a corresponding Jacobian 
matrix, which is rarely available in practice and is com¬ 
monly replaced with a forward difference (ED) approx¬ 
imation; see, e.g., [5]. Such approximate Newton-type 
optimization schemes utilize the ED approximation of 
the original nonlinear equation during every time step. 
An efficient variant of the approximate Newton-type op¬ 
timization can be performed by a Continuation NMFC 
(CNMFC) numerical method proposed by T. Ohtsuka 
in [8] , where each step of the algorithm requires solving 
a system of linear equations performed by the GMRES 
iterative method [10]. 

Our contributions presented below are two-fold. We 
describe an extension of CNMFC with a terminal con¬ 
straint, suitable to solve minimum-time optimal control 
problems, and with an optimization parameter. We in¬ 
vestigate preconditioning for GMRES in the context 
of the NMFC problems and using the MINRES iter¬ 
ation [9] instead of GMRES. MINRES provides overall 
faster implementation, compared to GMRES without 
restarts, of our approach in cases, where many iter¬ 
ations are required. Our numerical simulations show 
that the preconditioning can considerably improve the 
quality of controllers with marginal extra computational 
time, which can be reduced or eliminated by employing 
a parallel processing for the preconditioner setup. 

The rest of the paper is organized as follows. In Sec¬ 
tion 2, we formulate CNMFC of Ohtsuka, extended to 
having a terminal constraint and a parameter. Section 3 
describes the original algorithm of Ohtsuka, where the 
ED linear system is solved using GMRES, and then in¬ 
troduces MINRES as an alternative to GMRES, dis¬ 
cusses preconditioning for GMRES and MINRES, and 
suggests specific algorithms of constructing the precon¬ 
ditioner and using it to accelerate convergence of iter¬ 
ations. In Section 4, we give a detailed description of 
a test minimum-time optimal control problem, defin¬ 
ing a quickest arrival of the system to a given destina¬ 
tion, with inequality constraints on the system control, 
and its CNMFG formulation. Section 5 presents our re¬ 
sults of numerical experiments solving the test problem, 
demonstrating advantages of the proposed approaches. 



2 Finite horizon optimization by CNMPC 

As a specific example of a mathematical formalism of 
NMPC, we consider an extended version of the control 
problem considered by T. Ohtsuka [8] as follows, 

min J, 

u,p 

/•t+T 

J = (j){t+ T,x{t+ T),p) + J L{t',x{t'),u{t'),p)dt' 

subject to 

flT 

(2.2) C{t',x{t'),u{t'),p) =0, 


(2.3) ^{t + T,x{t + T),p) = 0. 

Here, x = x{t) denotes the vector of the state of the 
dynamic system, also serving as an initial state for the 
optimal control problem over the horizon. The vector 
u = u{t) is the control vector, serving as an input to 
control the system. The scalar function J describes 
a performance cost to be minimized, which includes a 
terminal cost (the first term in the sum) and a cost 
over the finite horizon (the second term in the sum). 
Equation (2.1) is the system dynamic model that may 
be nonlinear in x and/or u. Equation (2.2) describes the 
equality constraints for the state x and the control u. 
The horizon time length T may in principle also depend 
on t, e.g., for time-optimal control problems. In this 
case, the original problem can be converted into a fixed 
horizon problem by letting T{t) = 1 -tf, where </ is an 
additional parameter to be included in p and determined 
in MFC. Substituting t + rtf for the time t', we arrive 
at a problem with the normalized time scale r and fixed 
horizon Such a conversion is applied to the test 

problem in Section 4. 

Compared to [8], one extra constraint (2.3), de¬ 
scribed by the terminal constraint function tp, and an 
extra parameter vector p are being added to the prob¬ 
lem formulation, allowing one to extend CNMPC to a 
wide range of optimal control and design problems. 

The NMPC optimal control problem is solved by a 
variational approach. Its discrete counterpart is solved 
by the traditional Lagrange method of undetermined 
multipliers. We denote the costate vector by A and the 
Lagrange multiplier vector associated with the equality 
constraint (2.2) by p. The terminal constraint (2.3) is 
relaxed by introducing the Lagrange multiplier v. The 


so-called Hamiltonian function, as defined in control 
theory, is 

Hit, X, A, u, p, p) = L{t, X, u, p) 

-I- f{t,x,u,p) + p^C{t,x,u,p). 

To discretize the continuous formulation of the op¬ 
timal control problem stated above, we introduce a uni¬ 
form horizon time grid by dividing the horizon [t, t + T] 
into N time steps of size At and replace the time- 
continuous vector functions x{t) and u{t) by their in¬ 
dexed values Xi and Ui at the grid points. Thus, N is 
a number of artificial time steps for the optimal control 
problem over the horizon. The integral in the perfor¬ 
mance cost J over the time horizon is approximated by 
a simple quadrature rule. The time derivative of the 
state vector is approximated by the forward difference 
formula. Then the discretized optimal control problem 
appears as follows, 

min J, 

Ui ,p 

N-1 

J = (j){TN,XN,p) + ^ L{Ti,Xi,Ui,p)AT, 
i=0 

subject to 

Xi+i = Xi + f{Ti,Xi,Ui,p)AT, i = 0,l,...,N-1, 

C{Ti,Xi,Ui,p) =0, i = 0,1,..., A - 1, 

'iP{tn,xn:P) = 0 . 

We note that we have so far discretized the NMPC 
optimal control problem only in the horizon time. We 
will discretize the system time t later using the uniform 
time step size At, i.e. discretization in the horizon 
time may be different from the time discretization of 
the system. 

The necessary optimality conditions for the dis¬ 
cretized horizon problem are obtained using the discrete 
Lagrangian function 

N-l 

U) = (l){rN,XN:P) + ^ L{ri,Xi,Ui,p)AT 

i^O 

+ Ao[x(i) -Xo] 

N-l 

+ ^ >H-l-l[Xi -Xi^i f{Ti,Xi,Ui,p)AT] 
i=0 
N-l 

+ ^ fiJC{Ti,Xi,Ui,p)Ar -\-u'^'tl^{TN,XN,p), 

where X = [xi Ai]’^ and U = [ui pi v pY". Namely, 
the necessary optimality conditions coincide with the 
stationarity conditions 

^(A,C/)=0and^(A,[/) = 0. 




just obtained Xi and Ai, z = 0,1..., iV, as follows, 


For example, the derivative with respect to Ui, which is 
/dui = 0, yields the following equation: 


dL 

dui 


\T 

dui 

,dC 




+Mi -^{Ti,Xi,Ui,p)AT = 0. 


F[U,x,t] 

(ro, xo, A 1 , Mo, Mo, p) Ar 

^§^{Ti,Xi,X,+ i,Ui, tXi,p)AT 


Using the Hamiltonian function, it can be shortened to 


^^{tn-i,xn-i,Xn, mat-i, pn-i,p)At 


QH 

-—{Ti,Xi,Xi+l,Ui,Pi,p)AT = 0. 

UUi 


C{tq,xq,uq,p)At 


Taking the derivative with respect to pi, which is 
= 0, we obtain the following equation, which 
also involves the factor At, 


C{Ti,Xi,Ui,p)AT 


C (ttv -1, CCTV -1, MAT _ 1, p) At 


C{Ti,Xi,Ui,p)AT = 0. 

Now we proceed to the construction of a vector 
function F{U,x,t), which is used to formulate the full 
set of necessary optimality conditions. The vector 
function U = U{t) combines the control input m, the 
Lagrange multiplier p, the Lagrange multiplier v, and 
the parameter p, all in one vector, as follows, 

U{t) = [mJ, ..., ul_„p^, pI_„^^,pT- 

The vector argument x in the function F{U, x, t) denotes 
the current measured state vector, which serves as the 
initial vector xo in the following algorithm, defining an 
evaluation of F[U,x,t). 

1. Starting with the current measured state xg, com¬ 
pute Xi, i = 1,2 ... ,N, by the forward recursion 

Xi+i = Xi+ f{Ti,Xi,Ui,p)AT, i = Q,... ,N - 1. 

Then starting with the value 


tp{TN,XN,p) 

^{rN,XN:P) + ^{tn,XN:P)i^ 

dp {Xii Xi, Xi-\-i^Ui^ fJ,i,p)^T 

The optimality condition is the nonlinear equation 

(2.4) F[Uit),x{t),t]=0 

with respect to the unknown U(t), which needs to be 
solved numerically by a computer processor at each 
time step of NMPC in real time on the controller 
board. This is the most difficult and challenging part 
of implementation of NMPC. At the initial time t = to, 
we need to approximately solve (2.4) directly. 

Let us denote the step size of the system time 
discretization by At, assume that U{t — At) is already 
available at the time t, and set AU = U{t) — U{t — At). 
For a small scalar h > 0, which may be different from 
the system time step At and from the horizon time 
step At, we introduce the operator 

(2.5) a(y) = {F[U(t-At) + hV,x(t),t] 

— F[U{t — At),x{t),t])/h. 


Xn 


d(fy^ d'ljP' 

-^{tn,xn,p) + -^{tn,xn,p)v 


Then equation (2.4) is equivalent to the equation 
ha{AU/h) = b, where b = —F[U{t — At),x{t),t]. 


compute the costate Xi, i = N —1,... ,0, by the 
backward recursion 


Xi — Ai+i + 


dH^ 

dx 


{Ti,Xi, Xi+i,Ui, Pi,p)AT. 


Let us denote the j-th column of the mxm identity 
matrix by Cj, where m is the dimension of the vector U, 
and construct an m x m matrix A with the columns 
Acj, j = 1,... ,m, defined by the formula 

(2.6) Acj = a{ej). 


The matrix A approximates the symmetric Jacobian 
2. Calculate the vector function F[U,x,t], using the matrix Fu[U{t —At), x{t),t] so that a{V) = AV + 0{h). 







It is important to realize that the operator a(-) 
in (2.6) may be nonlinear. In particular, this explains 
why our algorithms of explicitly computing A for the 
purpose of a preconditioner setup may result in a non- 
symmetric matrix A. Numerical stability of computa¬ 
tions may be improved by enforcing the symmetry, by 
substituting {A + A’^)j2 for A. The deviation from the 
symmetry gets smaller with a sampling period h, which 
we are free to choose independently of At and At. 

A key limitation in the choice of h comes from the 
fact that the cancellation error starts picking up in the 
finite difference evaluation in the operator a(y) due to 
inexact arithmetic of the controller processor. This is 
an unavoidable side effect of using the finite difference 
approximation of the derivative. A recommended lower 
bound for the value of h can for example be 10“® in the 
double precision arithmetic, but the optimal value also 
depends on the function F[U,x,t]- 

Given the formulas for computing the vector func¬ 
tion F[U,x,t], nonlinear equation (2.4) must be solved 

at the points of the grid ti = to + iAt, z = 0,1,_ 

At the initial state xq = x{to), we hnd an approxi¬ 
mate solution Uo to the equation F[Uo,xoAo] = 0 by a 
suitable optimization procedure. The dimension of the 
vector u{t) is denoted by Uu- Since 

the first block entry of C/q, formed from the first 
elements of Uq, is taken as the control uq at the state xo- 
The next state xi = x{ti) is either measured by a sensor 
or computed by the formula xi = xq + Atf(to,xo,uo); 
cf. (2.1). Now we start the recursion as follows. 

At the time ti, where z > 0, we arrive with the state 
Xi and the vector Ui-i. The operator 

ai{V) = (F[Ui-i + hV, Xi, ti] - F[Ui-i,Xi, U]) jh, 

defined by (2.5), determines an m x m matrix Ai with 
the columns 

AiCj = ai{ej), j = l,...,m, 

as in (2.6). At the current time ti, our goal is to solve 
the following equation 

(2.7) hai{AUi/h) = 6^, where bi = -F[Ui-i,Xi,ti]. 

Then we set Ui = Ui-i + AUi and choose the first 
rzu components of Ui as the control Ui. The next state 
Xi+i = x{ti+i) either comes from a sensor, estimated, 
or computed by the formula Xi+i = Xi + Atf{ti, Xi, Ui). 

Having the basic setup of CNMPC now described, 
leading to equation (2.7), next we discuss numerical 
solution of (2.7). Let us highlight that equation (2.7) 
is never solved exactly in practice, thus, a choice of an 
algorithm may greatly affect not only the performance 
of the controller, but also the computed control as well. 


3 Algorithms 

A direct way to solve (2.7) approximately is generating 
the matrix Ai and then solving the system of linear 
equations AiAUi = bi by, e.g., the Gaussian elimination. 

Another way is solving (2.7) by a suitable Krylov 
subspace iteration, e.g., by GMRES [10] or MINRES [9] 
methods, where we do not need to generate the ma¬ 
trix Ai explicitly. Namely, we simply use the operator 
ai(y) instead of computing the matrix-vector product 
AiV, for arbitrary vectors V', cf., [5, 6]. In his seminal 
paper [8], T. Ohtsuka uses the GMRES iteration. 

A typical implementation of the preconditioned 
GMRES without restarts is given by Algorithm 1, where 
Tr denotes an action of a precontioner T on a vector r, 
as explained below. The unpreconditioned GMRES, as 
in [8], simply uses z = r. We denote by 
the submatrix of FI with the entries Hij such that 
H i < *2 and ji <j< j 2 - 


Algorithm 1 Preconditioned GMRES without restarts 

Input: a(z;), b, xq, fc^ax, T 
Output: Solution x of a{x) = b 
1 : r = b- a{xo), z = Tr, P = \\z\\ 2 , vi = z/P 
2: for fc = 1, . . . ,fcmax do 
3: r = a{vk), z = Tr 

4 : Hi,k,k = [vi,...,Vk]'^Z 

5 : Z = Z-[vi,...,Vk]Hi:k,k 

6: Hk+l.k = ll^lb 

7 : Vk+l = z/\\z \\2 

8: end for 

9: z/ = arg miny||i7i,fe„^^+i,i,fc„^^z/ - [/3,0,... ,0]^||2 
10 : x = xoA[vi,...,Vk^^Jiy 


We emphasize that the operator ap-) may be non¬ 
linear, but approximates the symmetric Jacobian ma¬ 
trix Fu \Ui-i,Xi, ti]. This implies a slight deviation from 
the symmetry property V 2 apVi) = {ai[V 2 ))^Vi for ar¬ 
bitrary vectors Vi and V 2 - We assume that the deviation 
is small and propose applying the MINRES iteration to 
solve equation (2.7). 

When the operator ai(-) is linear and symmetric, 
the projected (fcmax + 1) x k^ax matrix H, constructed 
by GMRES without preconditioning, is tridiagonal. 
The MINRES method is then a special variant of 
GMRES, which makes use of the tridiagonal structure. 
The table below, adopted from [3], gives a comparison 
of computational complexities of MINRES and GMRES 
without preconditioning for solution of a linear system 
Ax = b with a symmetric m x m matrix A in terms 
of memory storage required by working vectors in the 
solvers and the number of floating-point operations. By 
tp we denote the work needed for evaluating ai{V). 






Solver 

Storage 

Work per iteration 

MINRES 

7m 

tp + 9m 

GMRES 

(^max 2)777 

tp + {kmax + 3)m -f , ™ 

^max 


If the matrix Ai gets ill-conditioned, the conver¬ 
gence of GMRES or MINRES may stagnate. The con¬ 
vergence can be improved by preconditioning. A ma¬ 
trix Ti that approximates the matrix A~^ and such 
that computing the product Tir for an arbitrary vector 
r is relatively easy, is referred to as a preconditioner. 
The preconditioning for the system of linear equations 
Ax = b with the preconditioner T formally replaces 
the original system Ax = h with the equivalent pre¬ 
conditioned linear system TAx = Tb. If the condition 
number k{TA) = ||TA|| ||A“^r“^|| of the matrix TA is 
small, convergence of iterative solvers for the precondi¬ 
tioned system can be fast. However, the convergence 
of the preconditioned GMRES, in contrast to that of 
the preconditioned MINRES with a symmetric positive 
definite preconditioner, is not necessarily determined by 
the condition number k{TA). Results on convergence of 
GMRES in a nonlinear case can be found in [1]. 

When the approximate solution Xk^^^ computed by 
GMRES after fcmax iterations is not accurate enough, 
it is very common to restart GMRES with xq equal 
to Xk^^^ instead of increasing the maximum number of 
iterations fcmax- Practical implementations of GMRES 
perform restarts. Restarts allow to cap the GMRES 
memory use to fcmax + 2 vectors, but may significantly 
slow down the convergence. In our tests, we apply 
GMRES without restarts for simplicity of presentation. 

To setup the preconditioner, the matrix Ai is com¬ 
puted at some time ti and then its LU factorization 
Ai = LU is computed, where L is a lower- and U is 
an upper-triangular matrix. The product Tr is mathe¬ 
matically given by Tr = C/“^(L“^r), but is computed 
by back-substitution, which is much cheaper than the 
computation of the inverses of L and U. The same pre¬ 
conditioner T is used in a number of subsequent grid 
points starting from tj. The computation of the matrix 
Ai requires m evaluations ai[ej), see (2.6), that can be 
efficiently implemented in parallel. 

The symmetry of the preconditioner T can be used 
to reduce the memory storage and processor work; see, 
e.g., [2]. Eor example, the factorization T = LDL"’", see 
e.g. [4], instead of the LU factorization allows us using 
only half of memory. The anti-triangular factorization 
from [7] may also reduce both the memory requirements 
and work in preconditioning. 

MINRES requires symmetric positive definite pre¬ 
conditioners such as in [12]. In our MINRES simula¬ 
tions, although not reported in Section 5 in details, we 
use the preconditioned MINRES-QLP method from [3]. 


4 Test problem 

In this section, we formulate a test nonlinear prob¬ 
lem called TfG below for brevity, which describes the 
minimum-time motion from a state (a;o,2/o) to a state 
(x/, 2 //) with an inequality constrained control. 

The problem TfC has the following components: 


• State vector: x = 


X 

y 


. Input: u = 


u 

Ud 


• Parameter variables: p = [tf], where tf denotes the 
length of the evaluation horizon. 


• Dynamics: x = f{x,u,p) 


{Ax + B) cosw 
{Ax -|- B) sinu 


• Gonstraints: C{x,u,p) = [(u — Cu)^ + ~ = 0, 

i.e., the control u always stays within the band 

Cu - ru < u < Cu + Xu). 


Terminal constraints: ijj{x,p) = 


= 0 


X-Xf 

y-vf 

(the state should pass through the point {xf,yf) 
at t = tf) 


• Objective function to minimize: 

rt+tj 


nt + tf 

J = <f>{x,p)+ / L{x,u,p)dt', 


where 


cj){x,p) =tf, L{x,u,p) =-WdUd 

(the state should arrive at {xf,yf) in the shortest 
time; the function L serves to stabilize the slack 
variable Ud) 

• Gonstants: A = B = 1, xq = y^ = Q, t^ = Q, 
=yf = 1) Cu = 0.8, ru = 0.2, Wd = 0.005. 

The components of the corresponding discretized 
problem on the horizon are given below: 


• the scaled horizon time {t — To)/tf € [0,1] substi¬ 
tutes the original horizon time r € [to,to + tf]; 

• the discretized scaled horizon time is thus Ti = i Ar, 
where z = 0,1,..., A, and At = 1/N; 


the participating variables are the state 


Xi 


the 


costate 


Ai,i 

A2.i 


the control 


Ui 

'^di 


the Lagrange 


multipliers pi and 


1^1 

1^2 























• the state is governed by the model equation 

J Xi+i = Xi-\- At [p [Axi + B) cos m], 

1 2/i+i =yi + At [p {Axi + B) sin m], 

where z = 0,1,..., — 1; 

• the costate is determined by the backward recur¬ 
sion (Ai,Ar = Vi, A2,7V = V 2 ) 


Ai,i = Ai,i+i 

-I- At [pA(cosMiAi,i+i -I- sinMiA2.i-i-i)], 

A2,i = A2,i+1, 

where i = iV — 1, iV — 2, ..., 0; 

the equation F(U,xo,to) = 0, where 

U = [uo, Ud,0, ■ ■ ■ , UN-l,Ud,N-l, 

Pq, ... ,pN-i,’^ 1 , 1 ^ 2 ,p], 

has the following rows from the top to the bottom: 
Arp [{Axi + B) (- sinuiAi,i+i + cos UiA 2 .i-i-i) 

-f 2 {Ui Cu) Pi\ — 0 

Arp [2piUdi - Wd] = 0 


{ Arp [{ui - Cu)^ + u\ -rl] = 0 

J Xn — Xr = 0 
\ yN-yr = 0 


N-l 

At{ + S)(cositiAi,i+i + sinuiA 2 ,i+i) 

2—0 

-I- Pi[{ui - CuY + u^di - rl] - WdUdi} -1-1 = 0. 


Substituting ppi for pi, prior to differentiating the La- 
grangian, leads to alternative simpler and more numeri¬ 
cally stable, as observed in our tests, formulas, as follows 


At [p{Axi + B) (-sinMiAi,i-i-i + cosMiA 2 ,i+i) 

-f 2 (^Ui Pi] — 0 

At [2piUdi - Wdp] = 0 


{ At [{u, - Cuf + uli -rl] =0 

( Xn — Xr =Q 

\ yN-yr = 0 

{ JV-1 

At[ + B)(cosMiAi.i+i + sinitiA 2 .i+i) 

i=0 

- WdUdi] -1-1 = 0. 

We use the latter formulas in our numerical exper¬ 
iments described in the next section. 


5 Numerical results 

In our numerical experiments with the TfC problem the 
system of linear equations (2.7) is solved by the GMRES 
method. We have also tested MINRES, obtaining the 
controls similar to those with GMRES, reported here. 
The number of evaluations of a{V) in GMRES does not 
exceed an a priori chosen parameter denoted by /cmax, 
the error tolerance is tol = 10“®. The sampling time in 
the evaluation horizon is At = 0.1, the sampling time 
of the simulation is At = 0.02, and h = 10“®. 

The preconditioners are constructed as follows. At 
the time instances t = jtp, j = 0 , 1 ,..., with an a 
priori chosen time increment tp we calculate all entries 
of the matrix A by (2.6) and its LU factorization 
A = LU by Gaussian elimination with partial pivoting. 
The computed factors L and U are then used in the 
preconditioner as follows Tr = U~^{L~^r) for all 
sampling points t = iAt in the interval [jtp, (j -I- l)tp). 

The whole set of simulations reported here consists 
of the following four cases: 

1. no preconditioning, fcmax = 10; 

2. preconditioning with tp = 0.2 sec, fcmax = 1; 

3. preconditioning with tp = 0.4 sec, fcmax = 2; 

4. preconditioning with tp = 0.4 sec, fcmax = 10. 

The computed results are similar in all reported 

cases. Figure 1 displays the typical CNMPG control 
u, within the constant constraints, and the time to 
destination tj, both as functions of the system time in 
seconds, shown at the horizontal axis. Figure 2 shows a 
typical system trajectory in the x-y plane. 

Figures 3~6 show the value of ||F||, which we want to 
be vanished, and the GMRES residual (the left vertical 
axis) and the number of the actually performed GMRES 
iterations (the right vertical axis) at every system time 
step for all four cases, where the horizontal axis repre¬ 
sents the system time in seconds. Figure 3 corresponds 
to the GMRES iterations without preconditioning. Fig¬ 
ures 4-6 involve the preconditioner, recalculated with 
various frequencies, determined by the time increment 
tp, and for different fcmax ranging from 1 to 10. 

In Figure 3, the number of the actually performed 
GMRES iterations without preconditioning is always 
the maximum allowed in this test fcmax = 10. We use 
this test as a baseline for comparisons. 

We first point out a good behavior of the precon¬ 
ditioned GMRES even with fcmax = 1 and where the 
preconditioner is reconstructed once each tp = 0.2 sec, 
see Figure 4. This clearly demonstrates the fact that 
preconditioning reduces the number of evaluations of 
the vector function F(U,x,t). 


Residual Residual 
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— u 



°'^0 0.2 0.4 0.6 0.8 1 
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(1) NMPC control u and time to destination tf for TfC 
(reaches the target at t = 0.96) 



(3) GMRES without preconditioning, /cmax = 10 




(2) TfC trajectory by NPMC 



(4) GMRES with preconditionining, tp = 0.2 sec, fcmax = 1 





















































The effect of increasing the maximum number A:max 
of GMRES steps is seen by comparing Figures 4-6. 
Specihcally, in Figure 4, tp = 0.2 sec and fcmax = 1, 
compared to tp = 0.4 sec and fcmax = 2 in Figure 5, 
i.e., we can recompute the preconditioner twice less 
frequently at the cost of increasing fcmax from 1 to 2, 
and we observe a slightly better quality of the solution, 
as measured by the generally smaller values of ||F|| and 
the GMRES residual (the left vertical axis). 

In Figure 6, the preconditioner is recomputed as 
frequent as in Figure 5, but the largest allowed number 
of GMRES iterations is increased from fcmax = 2 to 
fcmax = 10. We observe in Figure 6 that GMRES often 
activates the default tolerance stopping criteria for the 
residual norm smaller than 10“^, before maxing out the 
allowed number of iterations fcmax- Overall, this leads to 
a generally much smaller residual in Figure 6 compared 
to that in Figure 5. However, the most decisive quantity 
11^11 behaves similar both in Figures 5 and 6, and the 
computed controls are so similar that the increase of 
fcmax from 2 to 10 may be unnecessary. 

Efficiency of preconditioning is illustrated by com¬ 
paring Figures 3 and 5, where the number of iterations is 
reduced five times giving similar/smaller values of HEH. 

In minimum-time optimal control problems, the 
length of the evaluation horizon gets smaller as the state 
(x, y) approaches the goal position. Near the goal posi¬ 
tion (1,1) the control has less capability (controllability) 
to direct the state towards the goal because of short time 
for control. This makes the equation F{U) = 0 more 
difficult for numerical solution, thus, ||E|| increases near 
the goal position, as seen in Figures 3-6. 

Conclusions 

Time-optimal problems are practically important, giv¬ 
ing optimal solutions for guidance, navigation and con¬ 
trol, which can be used for vehicles, trains, etc. Due 
to heavily nonlinear equations and highly coupled vari¬ 
ables, the time-optimal problems are difficult to solve 
numerically. We present an apparently first successful 
extension of CNMPC for real-time control of such prob¬ 
lems. Our numerical experiments demonstrate dramatic 
acceleration of convergence of iterations without sacri¬ 
ficing control quality, if proper preconditioning is used. 
The proposed concurrent construction of the precondi¬ 
tioner can be trivially efficiently implemented in parallel 
on controllers having multiple processing units, such as 
multi-core, graphics processing units, and modern field- 
programmable gate arrays. Replacing GMRES with the 
MINRES iterative solver may help reducing controller 
memory requirements and increasing the speed of con¬ 
vergence. Our algorithm, including the preconditioner 
setup implemented in parallel and the iterative solver, 


can significantly speed up the calculation of the control, 
compared to traditional sequential CNMPC algorithms, 
thus allowing to control system with faster dynamics. 
Our future work concerns analyzing MINRES, as a pos¬ 
sible replacement of GMRES, and developing efficient 
preconditioners, with faster on-line setup and applica¬ 
tion, within the framework of CNMPC. 
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